using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.QuadPack
{
    /// <summary>
    /// For further information view the QuadPack documentation.
    /// </summary>
    [Serializable]
    public static class DqextClass
    {
        /// <summary>
        /// For further information view the QuadPack documentation.
        /// </summary>
        /// <param name="n">The n parameter.</param>
        /// <param name="epstab">The epstab parameter.</param>
        /// <param name="abserr">The abserr parameter.</param>
        /// <param name="res3la">The res3la parameter.</param>
        /// <param name="nres">The nres parameter.</param>
        /// <returns>For further information view the QuadPack documentation.</returns>
        public static double Dqext(ref int n, double[] epstab, ref double abserr, double[] res3la, ref int nres)
        {
            double delta1;
            double delta2;
            double delta3;
            double epsinf;
            double error;
            double err1;
            double err2;
            double err3;
            double e0;
            double e1;
            double e1abs;
            double e2;
            double e3;
            double res;
            double result;
            double ss;
            double tol1;
            double tol2;
            double tol3;
            int NN;
            int i;
            int ib;
            int ib2;
            int ie;
            int indx;
            int k1;
            int k2;
            int k3;
            int limexp;
            int newelm;
            int num;

            (nres)++;
            NN = n;
            NN++; // make NN a FORTRAN array index 
            abserr = Double.MaxValue;
            result = epstab[n];
            if (NN < 3) // N < 3 
            {
                goto _100;
            }
            limexp = 50; // limexp = 50 
            epstab[n + 2] = epstab[n];
            newelm = (n) / 2; // (n-1)/2 
            epstab[n] = Double.MaxValue;
            num = NN;
            k1 = NN;
            for (i = 1; i <= newelm; i++)
            {
                k2 = k1 - 1;
                k3 = k1 - 2;
                res = epstab[k1 + 1];
                e0 = epstab[k3 - 1];
                e1 = epstab[k2 - 1];
                e2 = res;
                e1abs = Math.Abs(e1);
                delta2 = e2 - e1;
                err2 = Math.Abs(delta2);
                tol2 = (((Math.Abs(e2)) > (e1abs)) ? (Math.Abs(e2)) : (e1abs)) * Double.Epsilon;
                delta3 = e1 - e0;
                err3 = Math.Abs(delta3);
                tol3 = (((e1abs) > (Math.Abs(e0))) ? (e1abs) : (Math.Abs(e0))) * Double.Epsilon;
                if ((err2 > tol2) || (err3 > tol3))
                {
                    goto _10;
                }
                result = res;
                abserr = err2 + err3;
                goto _100;
                _10:
                e3 = epstab[k1 - 1];
                epstab[k1 - 1] = e1;
                delta1 = e1 - e3;
                err1 = Math.Abs(delta1);
                tol1 = (((e1abs) > (Math.Abs(e3))) ? (e1abs) : (Math.Abs(e3))) * Double.Epsilon;
                if ((err1 <= tol1) || (err2 <= tol2) || (err3 <= tol3))
                {
                    goto _20;
                }
                ss = 1.0 / delta1 + 1.0 / delta2 - 1.0 / delta3;
                epsinf = Math.Abs(ss * e1);
                if (epsinf > 1.0e-4)
                {
                    goto _30;
                }
                _20:
                NN = i + i - 1;
                goto _50;
                _30:
                res = e1 + 1.0 / ss;
                epstab[k1 - 1] = res;
                k1 -= 2;
                error = err2 + Math.Abs(res - e2) + err3;
                if (error > (abserr))
                {
                    goto _40;
                }
                abserr = error;
                result = res;
                _40:
                ;
            }
            _50:
            if (NN == limexp)
            {
                NN = 2 * (limexp / 2) - 1;
            }
            ib = 1; // ib = 1 
            if (((num / 2) * 2) == num) // ib = 2 
            {
                ib = 2;
            }
            ie = newelm + 1;
            for (i = 1; i <= ie; i++)
            {
                ib2 = ib + 2;
                epstab[ib - 1] = epstab[ib2 - 1];
                ib = ib2;
            }
            if (num == NN)
            {
                goto _80;
            }
            indx = num - NN + 1;
            for (i = 1; i <= NN; i++)
            {
                epstab[i - 1] = epstab[indx - 1];
                indx++;
            }
            _80:
            if (nres > 3) // nres >= 4 
            {
                goto _90;
            }
            res3la[(nres) - 1] = result;
            abserr = Double.MaxValue;
            goto _100;
            _90:
            abserr = Math.Abs(result - res3la[2]) + Math.Abs(result - res3la[1]) + Math.Abs(result - res3la[0]);
            res3la[0] = res3la[1];
            res3la[1] = res3la[2];
            res3la[2] = result;
            _100:
            abserr = (((abserr) > (5.0 * Double.Epsilon * Math.Abs(result)))
                          ? (abserr)
                          : (5.0 * Double.Epsilon * Math.Abs(result)));
            n = NN - 1;
            return result;
        }
    }
}